#!/usr/bin/perl
# 2013/05/08 Self ligation filter for ChIA-PET


use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;
use File::Basename;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -o [output filtered map file]\n";
}

my %opt;
getopts("i:o:", \%opt);
my $FILE_map = $opt{i};
my $FILE_out = $opt{o};
my ($name, $dir, $ext) = &fileparse($FILE_map, '\..*');
my $FILE_tmp = $dir . 'tmp_' . $name . '.txt';

my $COUNTER_self = 0;
my $COUNTER_other = 0;


#---------------------------------------
# keep self ligation candidate to rescue later
#---------------------------------------
my $fh_tmp = IO::File->new($FILE_tmp, 'w') or die "cannot write $FILE_tmp: $!";
my %data;

# id chr1 loc1 direction1 quality1 hind1 hind_loc1 chr2 loc2 direction2 quality2 hind2 hind_loc2 type(unique or repeat) number(alignable locations) locations
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;


	if(&checkSelfOrientation($chr1, $chr2, $loc1, $loc2, $direction1, $direction2) == 1){
		if(($chr1 eq $chr2) and (abs($loc1 - $loc2) < 10000)){
			$fh_tmp->print("$_\n");
		}else{
			$COUNTER_other++;
			$fh_out->print("$_\n");
		}
	}else{
		$fh_out->print("$_\n");
		$COUNTER_other++;
		if(($chr1 eq $chr2) and (abs($loc1 - $loc2) < 10200)){
			my $loc1_c = int($loc1 / 100) * 100;
			my $loc2_c = int($loc2 / 100) * 100;
			my ($key1, $key2);
			if($loc1_c < $loc2_c){
				($key1, $key2) = ("$chr1\t$loc1_c", "$chr2\t$loc2_c");
			}else{
				($key1, $key2) = ("$chr2\t$loc2_c", "$chr1\t$loc1_c");
			}
			$data{$key1}{$key2}++;
		}
	}
}
$fh_map->close();
$fh_tmp->close();


#---------------------------------------
# rescue from self ligation candidate
#---------------------------------------
my %addedNumber;
my $fh_rescue = IO::File->new($FILE_tmp) or die "cannot open $FILE_tmp: $!";
while($_ = $fh_rescue->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;
	my $loc1_c = int($loc1 / 100) * 100;
	my $loc2_c = int($loc2 / 100) * 100;

	my ($key1, $key2);
	if($loc1_c < $loc2_c){
		($key1, $key2) = ("$chr1\t$loc1_c", "$chr2\t$loc2_c");
	}else{
		($key1, $key2) = ("$chr2\t$loc2_c", "$chr1\t$loc1_c");
	}

	if(exists $data{$key1}{$key2}){
		my $Limit = $data{$key1}{$key2} / 3;
		if(exists $addedNumber{$key1}{$key2}){
			if($Limit - $addedNumber{$key1}{$key2} > 0.6){
				$COUNTER_other++;
				$fh_out->print("$_\n");
				$addedNumber{$key1}{$key2}++;
			}else{
				$COUNTER_self++;
			}
		}else{
			if($Limit > 0.6){
				$COUNTER_other++;
				$fh_out->print("$_\n");
				$addedNumber{$key1}{$key2} = 1;
			}
		}
	}else{
		$COUNTER_self++;
	}
}
$fh_rescue->close();
$fh_out->close();

unlink $FILE_tmp;

my $SUM = $COUNTER_other + $COUNTER_self;
print "Self ligation\n";
printf "Self ligation: %d (%.3f)\n", $COUNTER_self, $COUNTER_self / $SUM * 100;
printf "Non self ligation: %d (%.3f)\n\n", $COUNTER_other, $COUNTER_other / $SUM * 100;


sub checkSelfOrientation{
	my ($chr1, $chr2, $loc1, $loc2, $direction1, $direction2) = @_;
	if($chr1 ne $chr2){
		return 0;
	}
	if(($direction1 eq '+') and ($direction2 eq '+') and ($loc2 < $loc1)){
		return 1;
	}
	if(($direction1 eq '-') and ($direction2 eq '-') and ($loc1 < $loc2)){
		return 1;
	}
	return 0
}

